3. quality control
# protein remove outliners
tmp = pro.whole
pp = prcomp(t(tmp),cor=F)
outlinerids = c()
thetissues = unique(pro.whole.info$tissue_en)
for(i in 1:length(thetissues)){
idx = pro.whole.info$tissue_en == thetissues[i]
tmpinfo = pro.whole.info[idx,]
outx = is.outliner(pp$x[idx,1])
outy = is.outliner(pp$x[idx,2])
#outy = is.outliner(pp$x[idx,2],coef = 3)
if(sum(outx | outy) > 0){
outlinerids = c(outlinerids,rownames(tmpinfo)[outx | outy])
}
}
outlinerids
## [1] "X06080_Skin_of_back" "X11062_Pituitary" "X94356_Liver"
## [4] "X12092_Thyroid_gland" "X06080_Thyroid_gland" "X94356_Cecum"
## [7] "X16002_Adrenal_gland" "X92338_Adrenal_gland" "X12390_Fallopian_tube"
## [10] "X16068_Hypothalamus" "X06070_Pancreas" "X16086_Uterus"
vid = !is.element(colnames(pro.whole),outlinerids)
pro.whole = pro.whole[,vid]
pro.whole.info = pro.whole.info[vid,]
# met remove outliners
tmp = met.whole
pp = prcomp(t(tmp),cor=F)
outlinerids = c()
thetissues = unique(met.whole.info$tissue_en)
for(i in 1:length(thetissues)){
idx = met.whole.info$tissue_en == thetissues[i]
tmpinfo = met.whole.info[idx,]
outx = is.outliner(pp$x[idx,1])
outy = is.outliner(pp$x[idx,2])
#outy = is.outliner(pp$x[idx,2],coef = 3)
if(sum(outx | outy) > 0){
outlinerids = c(outlinerids,rownames(tmpinfo)[outx | outy])
}
}
outlinerids
## NULL
vid = !is.element(colnames(met.whole),outlinerids)
met.whole = met.whole[,vid]
met.whole.info = met.whole.info[vid,]
# mrna remove outliners
tmp = mrna.whole
pp = prcomp(t(tmp),cor=F)
outlinerids = c()
thetissues = unique(mrna.whole.info$tissue_en)
for(i in 1:length(thetissues)){
idx = mrna.whole.info$tissue_en == thetissues[i]
tmpinfo = mrna.whole.info[idx,]
outx = is.outliner(pp$x[idx,1])
outy = is.outliner(pp$x[idx,2])
#outy = is.outliner(pp$x[idx,2],coef = 3)
if(sum(outx | outy) > 0){
outlinerids = c(outlinerids,rownames(tmpinfo)[outx | outy])
}
}
outlinerids
## NULL
vid = !is.element(colnames(mrna.whole),outlinerids)
mrna.whole = mrna.whole[,vid]
mrna.whole.info = mrna.whole.info[vid,]
4. set colors
alltissues = names(pro.tissues)
tissue.systems = c('Integumentary','Endocrine','Brain','Respiratory','Digestive',
'Cardiovascular','Cardiovascular','Brain','Digestive','Endocrine',
'Cardiovascular','Muscle','Reproductive','Digestive','Brain',
'Immune','Renal','Endocrine','Digestive','Reproductive',
'Digestive','Brain','Muscle','Immune','Integumentary','Endocrine',
'Brain','Immune','Cardiovascular','Reproductive')
names(tissue.systems) = alltissues
tissue.color = pal_npg()(10)
names(tissue.color) = levels(factor(tissue.systems))
tissue.systems
## Skin_of_back Pituitary Frontal_pole
## "Integumentary" "Endocrine" "Brain"
## Lung Liver Arteria_cruralis
## "Respiratory" "Digestive" "Cardiovascular"
## Femoral_vein Hippocampus Ileocecum
## "Cardiovascular" "Brain" "Digestive"
## Thyroid_gland Arteria_carotis Muscle
## "Endocrine" "Cardiovascular" "Muscle"
## Ovary Cecum Superior_temporal_gyrus
## "Reproductive" "Digestive" "Brain"
## Spleen Kidney Adrenal_gland
## "Immune" "Renal" "Endocrine"
## Duodenum Fallopian_tube Stomach
## "Digestive" "Reproductive" "Digestive"
## Hypothalamus Heart Thymus
## "Brain" "Muscle" "Immune"
## Facial_skin Pancreas Supramarginal_gyrus
## "Integumentary" "Endocrine" "Brain"
## Adipose Aortic_arch Uterus
## "Immune" "Cardiovascular" "Reproductive"
tissue.color
## Brain Cardiovascular Digestive Endocrine Immune
## "#E64B35FF" "#4DBBD5FF" "#00A087FF" "#3C5488FF" "#F39B7FFF"
## Integumentary Muscle Renal Reproductive Respiratory
## "#8491B4FF" "#91D1C2FF" "#DC0000FF" "#7E6148FF" "#B09C85FF"
mypal = pal_aaas()(10)
mypal
## [1] "#3B4992FF" "#EE0000FF" "#008B45FF" "#631879FF" "#008280FF" "#BB0021FF"
## [7] "#5F559BFF" "#A20056FF" "#808180FF" "#1B1919FF"
show_col(mypal)
6 Figure 2: tissue aging DEGs and GO
# sub routine plot volcano
plot_Volcano <- function(res2, title){
tmpup = res2$Pvalue
tmpup[is.na(tmpup)] = 1
tmpup[res2$log2FC < 0] = 1
sortid_up = sort.int(-log10(tmpup),decreasing = T,index.return = T)$ix
tmpid.up = res2$ID[sortid_up[1:5]]
tmpdown = res2$Pvalue
tmpdown[is.na(tmpdown)] = 1
tmpdown[res2$log2FC > 0] = 1
sortid_down = sort.int(-log10(tmpdown),decreasing = T,index.return = T)$ix
tmpid.down = res2$ID[sortid_down[1:5]]
vid = is.element(res2$ID,c(tmpid.up,tmpid.down))
tlab = res2$ID
tlab[!vid] = NA
keyvals <- rep('gray50', nrow(res2))
names(keyvals) <- rep('NS', nrow(res2))
keyvals[which(res2$log2FC > 0.58 & res2$Pvalue < 0.05)] <- "Brown"
names(keyvals)[which(res2$log2FC > 0.58 & res2$Pvalue < 0.05)] <- 'High'
keyvals[which(res2$log2FC < -0.58 & res2$Pvalue < 0.05)] <- "darkblue"
names(keyvals)[which(res2$log2FC < -0.58 & res2$Pvalue < 0.05)] <- 'Low'
p = EnhancedVolcano(res2,
lab = tlab,
x = 'log2FC',
y = 'Pvalue',
caption = NULL,
title = title,
border = 'full',
titleLabSize = 12,
FCcutoff = 0.58,
cutoffLineWidth = F,
axisLabSize = 12,
subtitle = NULL,
cutoffLineCol = 'white',
gridlines.minor = F,
gridlines.major = F,
xlab = bquote(~Log[2]~ 'fold change'),
pCutoff = 0.05,
colCustom = keyvals,
colAlpha = 4/5,
legendPosition = 'none',
legendLabSize = 5,
legendIconSize = 3,
drawConnectors = TRUE,
widthConnectors = 0.5,
pointSize = -0.3*log10(res2$Pvalue),labSize = 3,
colConnectors = 'black')
return(p)
}
list_to_matrix <- function(DEproFC,alltissues){
DEproFC_matrix = list()
for(i in 1:length(alltissues)){
tmp = matrix(DEproFC[[i]],1,length(DEproFC[[i]]))
tmp = as.data.frame(tmp)
colnames(tmp) = names(DEproFC[[i]])
DEproFC_matrix[[i]] = tmp
}
DEproFC_matrix = t(as.matrix(rbindlist(DEproFC_matrix,fill = T)))
colnames(DEproFC_matrix) = names(DEproFC)
#vid = rowSums(is.na(DEproFC_matrix)) < ncol(DEproFC_matrix)/2
#DEproFC_matrix = DEproFC_matrix[vid,]
return(DEproFC_matrix)
}
met.class_enrichment <- function(mets,annote){
require(clusterProfiler)
vmet = intersect(mets,rownames(annote))
fluxgmt = data.frame(ont = annote$sub_class,
gene = rownames(annote),stringsAsFactors = F)
Recon3D <- enricher(gene = vmet,
TERM2GENE=fluxgmt,
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff =1,
qvalueCutoff = 1
)
Recon3Dout = Recon3D@result
Recon3Dout$DB = rep('HMDBclass',dim(Recon3Dout)[1])
n = dim(Recon3Dout)[2]
Recon3Dout = Recon3Dout[,c(n,1:(n-1))]
return(Recon3Dout)
}
6.1 protein
DEproFC = list()
DEproPvalue = list()
DEproAging = list()
DEproFC.develop = list()
DEproPvalue.develop = list()
for(i in 1:length(alltissues)){
#Mfuzz
thistissue = alltissues[i]
thispro = pro.tissues[[thistissue]]
thispro.header = pro.tissues.header[[thistissue]]
thispro = delete_dup_genes_forprotein(thispro,pro.tissues.header[[thistissue]])
thispro.info = pro.tissues.info[[thistissue]]
thisDEpro = DEGenes.simplified(thispro,catagory = thispro.info$stage == 4,
subset = thispro.info$stage == 4 | thispro.info$stage == 1)
thisDEpro.develop = DEGenes.simplified(thispro,catagory = thispro.info$stage == 2,
subset = thispro.info$stage == 2 | thispro.info$stage == 1)
forsort = thisDEpro$Pvalue
forsort[is.na(forsort)] = 1
idx = sort.int(forsort,decreasing = F,index.return = T)$ix
DEproAging[[i]] = thisDEpro[idx,-5]
DEproAging[[i]]$log2FC.devControl = thisDEpro.develop$log2FC
DEproAging[[i]]$Pvalue.devControl = thisDEpro.develop$Pvalue
DEproFC[[i]] = thisDEpro$log2FC
names(DEproFC[[i]]) = rownames(thisDEpro)
DEproPvalue[[i]] = thisDEpro$Pvalue
names(DEproPvalue[[i]]) = rownames(thisDEpro)
DEproFC.develop[[i]] = thisDEpro.develop$log2FC
names(DEproFC.develop[[i]]) = rownames(thisDEpro.develop)
DEproPvalue.develop[[i]] = thisDEpro.develop$Pvalue
names(DEproPvalue.develop[[i]]) = rownames(thisDEpro.develop)
}
names(DEproAging) = alltissues
names(DEproFC) = alltissues
names(DEproPvalue) = alltissues
names(DEproFC.develop) = alltissues
names(DEproPvalue.develop) = alltissues
proteinVoconoPlot = list()
names.DEproAging = names(DEproAging)
for(i in 1:length(DEproAging)){
res2 = DEproAging[[i]]
proteinVoconoPlot[[i]] = plot_Volcano(res2,names.DEproAging[i])
}
pdf(file = './out/20230217_aging/Figure2_DEG_GO_tissue/
Figure2A_DEpro_each_tissueV1.pdf',width = 3*5+1,height =3*6)
grid.arrange(arrangeGrob(grobs = proteinVoconoPlot,ncol = 5))
dev.off()
## png
## 2
#Data S2
openxlsx::write.xlsx(DEproAging, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S2_DE_tissue_Aging_pro.xlsx")
# reducce size Data S2
tpath = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S2_DE_tissue_Aging_pro.xlsx"
sheetNames = openxlsx::getSheetNames(tpath)
xx = list()
for(i in 1:length(sheetNames)){
tmp = openxlsx::read.xlsx(tpath,sheet = sheetNames[i])
tmp[c(2,3,4,7,8)] = signif(tmp[c(2,3,4,7,8)],3)
tmp = tmp[,-c(5,6)]
xx[[i]] = tmp
}
names(xx) = sheetNames
openxlsx::write.xlsx(xx, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/reduce_Data S2_DE_tissue_Aging_pro.xlsx")
# num up and down proteins
DEproFC_matrix = list_to_matrix(DEproFC,alltissues)
DEproPvalue_matrix = list_to_matrix(DEproPvalue,alltissues)
DEproFC_matrix.develop = list_to_matrix(DEproFC.develop,alltissues)
DEproPvalue_matrix.develop = list_to_matrix(DEproPvalue.develop,alltissues)
Aging_pro_sigup_matrix = (DEproFC_matrix > 0.58 & DEproPvalue_matrix < 0.05
& DEproPvalue_matrix.develop > 0.05) +0
Aging_pro_sigdown_matrix = -((DEproFC_matrix < -0.58 & DEproPvalue_matrix < 0.05
& DEproPvalue_matrix.develop > 0.05) +0)
Aging_pro_sigall_matrix = Aging_pro_sigup_matrix + Aging_pro_sigdown_matrix
Aging_pro_updown = data.frame(stringsAsFactors = F,
num.up = colSums(Aging_pro_sigup_matrix,na.rm =T)/
colSums(!is.na(Aging_pro_sigup_matrix)),
num.down = colSums(Aging_pro_sigdown_matrix,na.rm =T)/
colSums(!is.na(Aging_pro_sigdown_matrix)),
num.all = colSums(abs(Aging_pro_sigall_matrix),na.rm =T)/
colSums(!is.na(Aging_pro_sigall_matrix)),
tissues = colnames(Aging_pro_sigup_matrix),
tissue_systems = tissue.systems)
rownames(Aging_pro_updown) = Aging_pro_updown$tissues
Aging_pro_updown
## num.up num.down num.all
## Skin_of_back 0.036514823 -0.02422270 0.06081448
## Pituitary 0.024935277 -0.02820548 0.05316973
## Frontal_pole 0.022538363 -0.01742606 0.03996803
## Lung 0.033295982 -0.01888302 0.05220176
## Liver 0.034220532 -0.01801802 0.05225225
## Arteria_cruralis 0.025117739 -0.04575335 0.07101526
## Femoral_vein 0.041103448 -0.03337010 0.07450331
## Hippocampus 0.013891271 -0.02023320 0.03413379
## Ileocecum 0.021394879 -0.04090909 0.06238408
## Thyroid_gland 0.037166086 -0.01646011 0.05364059
## Arteria_carotis 0.024580336 -0.02637890 0.05096942
## Muscle 0.062205062 -0.06260720 0.12483912
## Ovary 0.070639717 -0.08727348 0.15794957
## Cecum 0.040690691 -0.01725431 0.05795796
## Superior_temporal_gyrus 0.015867713 -0.01687270 0.03274307
## Spleen 0.019508449 -0.03287250 0.05238900
## Kidney 0.018391573 -0.01220532 0.03060201
## Adrenal_gland 0.042047532 -0.01391427 0.05604055
## Duodenum 0.030563661 -0.06327753 0.09387003
## Fallopian_tube 0.029572113 -0.04311911 0.07273930
## Stomach 0.027610674 -0.05751735 0.08519833
## Hypothalamus 0.011605416 -0.00902498 0.02063185
## Heart 0.043579314 -0.02086957 0.06451613
## Thymus 0.075391850 -0.32503133 0.40075259
## Facial_skin 0.009045226 -0.01404917 0.02313301
## Pancreas 0.011964948 -0.01736638 0.02933738
## Supramarginal_gyrus 0.017594835 -0.01743904 0.03504522
## Adipose 0.008937121 -0.08154781 0.09053103
## Aortic_arch 0.066784870 -0.07934633 0.14626454
## Uterus 0.030058007 -0.04924376 0.07935949
## tissues tissue_systems
## Skin_of_back Skin_of_back Integumentary
## Pituitary Pituitary Endocrine
## Frontal_pole Frontal_pole Brain
## Lung Lung Respiratory
## Liver Liver Digestive
## Arteria_cruralis Arteria_cruralis Cardiovascular
## Femoral_vein Femoral_vein Cardiovascular
## Hippocampus Hippocampus Brain
## Ileocecum Ileocecum Digestive
## Thyroid_gland Thyroid_gland Endocrine
## Arteria_carotis Arteria_carotis Cardiovascular
## Muscle Muscle Muscle
## Ovary Ovary Reproductive
## Cecum Cecum Digestive
## Superior_temporal_gyrus Superior_temporal_gyrus Brain
## Spleen Spleen Immune
## Kidney Kidney Renal
## Adrenal_gland Adrenal_gland Endocrine
## Duodenum Duodenum Digestive
## Fallopian_tube Fallopian_tube Reproductive
## Stomach Stomach Digestive
## Hypothalamus Hypothalamus Brain
## Heart Heart Muscle
## Thymus Thymus Immune
## Facial_skin Facial_skin Integumentary
## Pancreas Pancreas Endocrine
## Supramarginal_gyrus Supramarginal_gyrus Brain
## Adipose Adipose Immune
## Aortic_arch Aortic_arch Cardiovascular
## Uterus Uterus Reproductive
mean(Aging_pro_updown$num.all)
## [1] 0.07529843
6.2 mRNA
DEmrnaFC = list()
DEmrnaPvalue = list()
DEmrnaAging = list()
DEmrnaFC.develop = list()
DEmrnaPvalue.develop = list()
for(i in 1:length(alltissues)){
#Mfuzz
thistissue = alltissues[i]
thismrna = mrna.tissues[[thistissue]]
#thismrna.header = mrna.tissues.header[[thistissue]]
thismrna.info = mrna.tissues.info[[thistissue]]
if (sum(thismrna.info$stage == 1) ==1){
thismrna = cbind(thismrna[,thismrna.info$stage == 1],thismrna)
thismrna.info = rbind(thismrna.info[thismrna.info$stage == 1,],thismrna.info)
}
if (sum(thismrna.info$stage == 4) ==1){
thismrna = cbind(thismrna,thismrna[,thismrna.info$stage == 4])
thismrna.info = rbind(thismrna.info,thismrna.info[thismrna.info$stage == 4,])
}
thisDEmrna = DEGenes.simplified(thismrna,catagory = thismrna.info$stage == 4,
subset = thismrna.info$stage == 4 | thismrna.info$stage == 1)
thisDEmrna.develop = DEGenes.simplified(thismrna,catagory = thismrna.info$stage == 2,
subset = thismrna.info$stage == 2 | thismrna.info$stage == 1)
idx = sort.int(thisDEmrna$Pvalue,decreasing = F,index.return = T)$ix
DEmrnaAging[[i]] = thisDEmrna[idx,-5]
DEmrnaAging[[i]]$log2FC.devControl = thisDEmrna.develop$log2FC
DEmrnaAging[[i]]$Pvalue.devControl = thisDEmrna.develop$Pvalue
DEmrnaFC[[i]] = thisDEmrna$log2FC
names(DEmrnaFC[[i]]) = rownames(thisDEmrna)
DEmrnaPvalue[[i]] = thisDEmrna$Pvalue
names(DEmrnaPvalue[[i]]) = rownames(thisDEmrna)
DEmrnaFC.develop[[i]] = thisDEmrna.develop$log2FC
names(DEmrnaFC.develop[[i]]) = rownames(thisDEmrna.develop)
DEmrnaPvalue.develop[[i]] = thisDEmrna.develop$Pvalue
names(DEmrnaPvalue.develop[[i]]) = rownames(thisDEmrna.develop)
}
names(DEmrnaFC) = alltissues
names(DEmrnaPvalue) = alltissues
names(DEmrnaAging) = alltissues
names(DEmrnaFC.develop) = alltissues
names(DEmrnaPvalue.develop) = alltissues
DEmrnaFC_matrix = list_to_matrix(DEmrnaFC,alltissues)
DEmrnaPvalue_matrix = list_to_matrix(DEmrnaPvalue,alltissues)
DEmrnaFC_matrix.develop = list_to_matrix(DEmrnaFC.develop,alltissues)
DEmrnaPvalue_matrix.develop = list_to_matrix(DEmrnaPvalue.develop,alltissues)
#data S1
openxlsx::write.xlsx(DEmrnaAging, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S1_DE_tissue_Aging_mRNA.xlsx")
# reducce size Data S1
tpath = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S1_DE_tissue_Aging_mRNA.xlsx"
sheetNames = openxlsx::getSheetNames(tpath)
xx = list()
for(i in 1:length(sheetNames)){
tmp = openxlsx::read.xlsx(tpath,sheet = sheetNames[i])
tmp[c(2,3,4,7,8)] = signif(tmp[c(2,3,4,7,8)],3)
tmp = tmp[,-c(5,6)]
xx[[i]] = tmp
}
names(xx) = sheetNames
openxlsx::write.xlsx(xx, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/reduce_Data S1_DE_tissue_Aging_mRNA.xlsx")
Aging_mrna_sigup_matrix = (DEmrnaFC_matrix > 0.58 & DEmrnaPvalue_matrix < 0.05 &
DEmrnaPvalue_matrix.develop > 0.05) +0
Aging_mrna_sigdown_matrix = -((DEmrnaFC_matrix < -0.58 & DEmrnaPvalue_matrix < 0.05 &
DEmrnaPvalue_matrix.develop > 0.05) +0)
Aging_mrna_sigall_matrix = Aging_mrna_sigup_matrix + Aging_mrna_sigdown_matrix
Aging_mrna_updown = data.frame(stringsAsFactors = F,
num.up = colSums(Aging_mrna_sigup_matrix,na.rm =T)/
colSums(!is.na(Aging_mrna_sigup_matrix)),
num.down = colSums(Aging_mrna_sigdown_matrix,na.rm =T)/
colSums(!is.na(Aging_mrna_sigdown_matrix)),
num.all = colSums(abs(Aging_mrna_sigall_matrix),na.rm =T)/
colSums(!is.na(Aging_mrna_sigall_matrix)),
tissues = colnames(Aging_mrna_sigup_matrix),
tissue_systems = tissue.systems)
rownames(Aging_mrna_updown) = Aging_mrna_updown$tissues
Aging_mrna_updown
## num.up num.down num.all
## Skin_of_back 0.014877204 -0.013853904 0.028731108
## Pituitary 0.016756193 -0.020444159 0.037200353
## Frontal_pole 0.017589731 -0.020362887 0.037952619
## Lung 0.008547009 -0.015562006 0.024109015
## Liver 0.015846995 -0.011202186 0.027049180
## Arteria_cruralis 0.016141630 -0.013364575 0.029506205
## Femoral_vein 0.012257843 -0.012772158 0.025030002
## Hippocampus 0.020030204 -0.019791749 0.039821954
## Ileocecum 0.054235073 -0.028506085 0.082741158
## Thyroid_gland 0.026949541 -0.013761468 0.040711009
## Arteria_carotis 0.021807382 -0.029953331 0.051760713
## Muscle 0.014691302 -0.028301037 0.042992339
## Ovary 0.025969730 -0.035232818 0.061202547
## Cecum 0.050902852 -0.040853011 0.091755862
## Superior_temporal_gyrus 0.015590200 -0.023783010 0.039373210
## Spleen 0.029801597 -0.034164814 0.063966411
## Kidney 0.020382901 -0.013863674 0.034246575
## Adrenal_gland 0.017783149 -0.026588398 0.044371547
## Duodenum 0.119320039 -0.126266754 0.245586793
## Fallopian_tube 0.049622774 -0.034300381 0.083923155
## Stomach 0.009763363 -0.008439517 0.018202879
## Hypothalamus 0.001815598 -0.003473319 0.005288917
## Heart 0.009539101 -0.011616925 0.021156026
## Thymus 0.006846780 -0.001572909 0.008419689
## Facial_skin 0.039662514 -0.040293329 0.079955843
## Pancreas 0.024070432 -0.024445069 0.048515501
## Supramarginal_gyrus 0.014104710 -0.022870348 0.036975058
## Adipose 0.014931669 -0.013497554 0.028429222
## Aortic_arch 0.029785810 -0.027024766 0.056810576
## Uterus 0.020403106 -0.025028911 0.045432017
## tissues tissue_systems
## Skin_of_back Skin_of_back Integumentary
## Pituitary Pituitary Endocrine
## Frontal_pole Frontal_pole Brain
## Lung Lung Respiratory
## Liver Liver Digestive
## Arteria_cruralis Arteria_cruralis Cardiovascular
## Femoral_vein Femoral_vein Cardiovascular
## Hippocampus Hippocampus Brain
## Ileocecum Ileocecum Digestive
## Thyroid_gland Thyroid_gland Endocrine
## Arteria_carotis Arteria_carotis Cardiovascular
## Muscle Muscle Muscle
## Ovary Ovary Reproductive
## Cecum Cecum Digestive
## Superior_temporal_gyrus Superior_temporal_gyrus Brain
## Spleen Spleen Immune
## Kidney Kidney Renal
## Adrenal_gland Adrenal_gland Endocrine
## Duodenum Duodenum Digestive
## Fallopian_tube Fallopian_tube Reproductive
## Stomach Stomach Digestive
## Hypothalamus Hypothalamus Brain
## Heart Heart Muscle
## Thymus Thymus Immune
## Facial_skin Facial_skin Integumentary
## Pancreas Pancreas Endocrine
## Supramarginal_gyrus Supramarginal_gyrus Brain
## Adipose Adipose Immune
## Aortic_arch Aortic_arch Cardiovascular
## Uterus Uterus Reproductive
mean(Aging_mrna_updown$num.all)
## [1] 0.04937392
mrnaVoconoPlot = list()
names.DEmrnaAging = names(DEmrnaAging)
for(i in 1:length(DEmrnaAging)){
res2 = DEmrnaAging[[i]]
mrnaVoconoPlot[[i]] = plot_Volcano(res2,names.DEmrnaAging[i])
}
#this plot is large, plot to file
pdf(file = './out/20230217_aging/Figure2_DEG_GO_tissue/FigureS2B_DEmrna_each_tissueV1.pdf',
width = 3*5+1,height =3*6)
grid.arrange(arrangeGrob(grobs = mrnaVoconoPlot,ncol = 5))
dev.off()
## png
## 2
6.3 metabolism
DEmetFC = list()
DEmetPvalue = list()
DEmetAging = list()
DEmetFC.develop = list()
DEmetPvalue.develop = list()
for(i in 1:length(alltissues)){
#Mfuzz
thistissue = alltissues[i]
thismet = met.tissues[[thistissue]]
thismet.header = met.tissues.header[[thistissue]]
thismet.info = met.tissues.info[[thistissue]]
if (sum(thismet.info$stage == 1) ==1){
thismet = cbind(thismet[,thismet.info$stage == 1],thismet)
thismet.info = rbind(thismet.info[thismet.info$stage == 1,],thismet.info)
}
if (sum(thismet.info$stage == 4) ==1){
thismet = cbind(thismet,thismet[,thismet.info$stage == 4])
thismet.info = rbind(thismet.info,thismet.info[thismet.info$stage == 4,])
}
thisDEmet = DEGenes.simplified(thismet,catagory = thismet.info$stage == 4,
subset = thismet.info$stage == 4 | thismet.info$stage == 1)
thisDEmet.develop = DEGenes.simplified(thismet,catagory = thismet.info$stage == 2,
subset = thismet.info$stage == 2 | thismet.info$stage == 1)
forsort = thisDEmet$Pvalue
forsort[is.na(forsort)] = 1
idx = sort.int(forsort,decreasing = F,index.return = T)$ix
DEmetAging[[i]] = thisDEmet[idx,-5]
DEmetAging[[i]]$log2FC.devControl = thisDEmet.develop$log2FC
DEmetAging[[i]]$Pvalue.devControl = thisDEmet.develop$Pvalue
DEmetFC[[i]] = thisDEmet$log2FC
names(DEmetFC[[i]]) = rownames(thisDEmet)
DEmetPvalue[[i]] = thisDEmet$Pvalue
names(DEmetPvalue[[i]]) = rownames(thisDEmet)
DEmetFC.develop[[i]] = thisDEmet.develop$log2FC
names(DEmetFC.develop[[i]]) = rownames(thisDEmet.develop)
DEmetPvalue.develop[[i]] = thisDEmet.develop$Pvalue
names(DEmetPvalue.develop[[i]]) = rownames(thisDEmet.develop)
}
names(DEmetAging) = alltissues
names(DEmetFC) = alltissues
names(DEmetPvalue) = alltissues
DEmetFC_matrix = list_to_matrix(DEmetFC,alltissues)
DEmetPvalue_matrix = list_to_matrix(DEmetPvalue,alltissues)
DEmetFC_matrix.develop = list_to_matrix(DEmetFC.develop,alltissues)
DEmetPvalue_matrix.develop = list_to_matrix(DEmetPvalue.develop,alltissues)
#write Data S3
openxlsx::write.xlsx(DEmetAging, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S3_DE_tissue_Aging_metabolite.xlsx")
# reducce size Data S2
tpath = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S3_DE_tissue_Aging_metabolite.xlsx"
sheetNames = openxlsx::getSheetNames(tpath)
xx = list()
for(i in 1:length(sheetNames)){
tmp = openxlsx::read.xlsx(tpath,sheet = sheetNames[i])
tmp[c(2,3,4,7,8)] = signif(tmp[c(2,3,4,7,8)],3)
tmp = tmp[,-c(5,6)]
xx[[i]] = tmp
}
names(xx) = sheetNames
openxlsx::write.xlsx(xx, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/reduce_Data S3_DE_tissue_Aging_metabolite.xlsx")
# num DEmets
Aging_met_sigup_matrix = (DEmetFC_matrix > 0.58 & DEmetPvalue_matrix < 0.05 &
DEmetPvalue_matrix.develop > 0.05) + 0
Aging_met_sigdown_matrix = -((DEmetFC_matrix < -0.58 & DEmetPvalue_matrix < 0.05 &
DEmetPvalue_matrix.develop > 0.05) + 0)
Aging_met_sigall_matrix = Aging_met_sigup_matrix + Aging_met_sigdown_matrix
Aging_met_updown = data.frame(stringsAsFactors = F,
num.up = colSums(Aging_met_sigup_matrix,na.rm =T)/
colSums(!is.na(Aging_met_sigup_matrix)),
num.down = colSums(Aging_met_sigdown_matrix,na.rm =T)/
colSums(!is.na(Aging_met_sigdown_matrix)),
num.all = colSums(abs(Aging_met_sigall_matrix),na.rm =T)/
colSums(!is.na(Aging_met_sigall_matrix)),
tissues = colnames(Aging_met_sigup_matrix),
tissue_systems = tissue.systems)
rownames(Aging_met_updown) = Aging_met_updown$tissues
Aging_met_updown
## num.up num.down num.all
## Skin_of_back 0.036030341 -0.027180784 0.06321113
## Pituitary 0.127226463 -0.078032231 0.20525869
## Frontal_pole 0.013303769 -0.047302291 0.06060606
## Lung 0.015033073 -0.021046302 0.03607937
## Liver 0.006362059 -0.029496819 0.03585888
## Arteria_cruralis 0.020116807 -0.020765737 0.04088254
## Femoral_vein 0.000000000 -0.011150235 0.01115023
## Hippocampus 0.014792899 -0.015532544 0.03032544
## Ileocecum 0.012651265 -0.058305831 0.07095710
## Thyroid_gland 0.006051437 -0.084720121 0.09077156
## Arteria_carotis 0.048211509 -0.010886470 0.05909798
## Muscle 0.011284722 -0.015625000 0.02690972
## Ovary 0.040852575 -0.080521018 0.12137359
## Cecum 0.147406266 -0.083204931 0.23061120
## Superior_temporal_gyrus 0.065759637 -0.077853364 0.14361300
## Spleen 0.053291536 -0.035736677 0.08902821
## Kidney 0.009730967 -0.062392673 0.07212364
## Adrenal_gland 0.021821632 -0.013757116 0.03557875
## Duodenum 0.034172662 -0.050959233 0.08513189
## Fallopian_tube 0.006607930 -0.042951542 0.04955947
## Stomach 0.006823821 -0.047766749 0.05459057
## Hypothalamus 0.030780781 -0.024024024 0.05480480
## Heart 0.028863500 -0.023451594 0.05231509
## Thymus 0.023887588 -0.136299766 0.16018735
## Facial_skin 0.002056555 -0.011825193 0.01388175
## Pancreas 0.026440678 -0.016949153 0.04338983
## Supramarginal_gyrus 0.126301180 -0.045801527 0.17210271
## Adipose 0.002421308 -0.007263923 0.00968523
## Aortic_arch 0.015799868 -0.048716261 0.06451613
## Uterus 0.031384615 -0.035692308 0.06707692
## tissues tissue_systems
## Skin_of_back Skin_of_back Integumentary
## Pituitary Pituitary Endocrine
## Frontal_pole Frontal_pole Brain
## Lung Lung Respiratory
## Liver Liver Digestive
## Arteria_cruralis Arteria_cruralis Cardiovascular
## Femoral_vein Femoral_vein Cardiovascular
## Hippocampus Hippocampus Brain
## Ileocecum Ileocecum Digestive
## Thyroid_gland Thyroid_gland Endocrine
## Arteria_carotis Arteria_carotis Cardiovascular
## Muscle Muscle Muscle
## Ovary Ovary Reproductive
## Cecum Cecum Digestive
## Superior_temporal_gyrus Superior_temporal_gyrus Brain
## Spleen Spleen Immune
## Kidney Kidney Renal
## Adrenal_gland Adrenal_gland Endocrine
## Duodenum Duodenum Digestive
## Fallopian_tube Fallopian_tube Reproductive
## Stomach Stomach Digestive
## Hypothalamus Hypothalamus Brain
## Heart Heart Muscle
## Thymus Thymus Immune
## Facial_skin Facial_skin Integumentary
## Pancreas Pancreas Endocrine
## Supramarginal_gyrus Supramarginal_gyrus Brain
## Adipose Adipose Immune
## Aortic_arch Aortic_arch Cardiovascular
## Uterus Uterus Reproductive
mean(Aging_met_updown$num.all)
## [1] 0.07502263
metVoconoPlot = list()
names.DEmetAging = names(DEmetAging)
for(i in 1:length(DEmetAging)){
res2 = DEmetAging[[i]]
metVoconoPlot[[i]] = plot_Volcano(res2,names.DEmetAging[i])
}
pdf(file = './out/20230217_aging/Figure2_DEG_GO_tissue/
FigureS2C_DEmet_each_tissueV1.pdf',width = 3*5+1,height =3*6)
grid.arrange(arrangeGrob(grobs = metVoconoPlot,ncol = 5))
dev.off()
## png
## 2
6.4 Figure 2A perc changed mols each tissue
fmt_dcimals <- function(decimals=0){
function(x) format(x,nsmall = decimals,scientific = FALSE)
}
plotAgingNum = list()
idx = sort.int(Aging_pro_updown$num.all,decreasing = F,index.return = T)$ix
Aging_pro_updown.v = Aging_pro_updown[idx,]
idx = sort.int(Aging_pro_updown.v$tissue_systems,decreasing = F,index.return = T)$ix
Aging_pro_updown.v = Aging_pro_updown.v[idx,]
tissueindex = rownames(Aging_pro_updown.v)
#metabolites
Aging_met_updown.v = Aging_met_updown[tissueindex,]
p1 = ggplot(Aging_met_updown.v,aes(x = factor(tissues,level = tissues),y = num.up,
color = tissue_systems,fill = tissue_systems))+
geom_bar(stat="identity",alpha = 0.8)
plotAgingNum[[1]] = p1+ geom_bar(stat="identity",aes(y = num.down),alpha = 0.6)+
geom_hline(yintercept = 0)+
theme_classic()+lghplot.addtheme(legend.position = 'none',hjust = 1,size = 13.5)+
theme(axis.text.y = element_text(size = 9.5, face = "bold", color = "black"))+
scale_color_aaas(alpha = 0.6)+
scale_fill_aaas(alpha = 0.6)+theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),axis.line.x =element_blank())+
scale_y_continuous(labels = fmt_dcimals(2))+
ylab('Metabolites')+ xlab('')
#protein
Aging_pro_updown.v = Aging_pro_updown[tissueindex,]
p1 = ggplot(Aging_pro_updown.v,aes(x = factor(tissues,level = tissues),y = num.up,
color = tissue_systems,fill = tissue_systems))+
geom_bar(stat="identity",alpha = 0.8)
plotAgingNum[[2]] = p1+ geom_bar(stat="identity",aes(y = num.down),alpha = 0.6)+
geom_hline(yintercept = 0)+
theme_classic()+lghplot.addtheme(legend.position = 'none',hjust = 1,size = 13.5)+
theme(axis.text.y = element_text(size = 9.5, face = "bold", color = "black"))+
scale_color_aaas(alpha = 0.6)+
scale_fill_aaas(alpha = 0.6)+
scale_y_continuous(labels = fmt_dcimals(2))+
ylab('Proteins')+ xlab('')
#mrna
Aging_mrna_updown.v = Aging_mrna_updown[tissueindex,]
p1 = ggplot(Aging_mrna_updown.v,aes(x = factor(tissues,level = tissues),y = num.up,
color = tissue_systems,fill = tissue_systems))+
geom_bar(stat="identity",alpha = 0.8)
plotAgingNum[[3]] = p1+ geom_bar(stat="identity",aes(y = num.down),alpha = 0.6)+
geom_hline(yintercept = 0)+
theme_classic()+lghplot.addtheme(legend.position = 'none',hjust = 1,size = 13.5)+
theme(axis.text.y = element_text(size = 9.5, face = "bold", color = "black"))+
scale_color_aaas(alpha = 0.6)+
scale_fill_aaas(alpha = 0.6)+scale_y_continuous(labels = fmt_dcimals(2))+
ylab('mRNAs')+ xlab('')
#pdf(file = './out/20230217_aging/Figure2_DEG_GO_tissue/
# Figure2A_number_of_Aging_moleculars_prometv1_aaas.pdf',width = 12,height =5)
grid.arrange(arrangeGrob(grobs = plotAgingNum[1:2],
top = 'Overall molecular changes with aging',
ncol = 1,heights = c(1.7,3.3)))
#dev.off()
6.5 Figure 2B GO enrichment
#write for metascape up
inputGOenrich = data.frame(tissues = names(DEproAging),stringsAsFactors = F,
genes = rep('c',length(DEproAging)))
#colnames(inputGOenrich) = c('tissue','GeneSymbol')
for(i in 1:length(DEproAging)){
tmp = DEproAging[[i]]
tgenes = unique(rownames(tmp)[tmp$Pvalue < 0.05 & tmp$log2FC > 0.58 &
tmp$Pvalue.devControl > 0.05])
tgenes = tgenes[!is.na(tgenes)]
tgenes = tgenes[tgenes != '']
tgenes = paste0(tgenes,collapse =',')
inputGOenrich$genes[i] = tgenes
}
writetxt(inputGOenrich,'./out/20230217_aging/Figure2_DEG_GO_tissue/protein/_inputGOenrich_age_DEpro_up.txt',row.names = F,col.names = F)
#write for metascape down
inputGOenrich = data.frame(tissues = names(DEproAging),stringsAsFactors = F,
genes = rep('c',length(DEproAging)))
#colnames(inputGOenrich) = c('tissue','GeneSymbol')
for(i in 1:length(DEproAging)){
tmp = DEproAging[[i]]
tgenes = unique(rownames(tmp)[tmp$Pvalue < 0.05 & tmp$log2FC < -0.58 &
tmp$Pvalue.devControl > 0.05])
tgenes = tgenes[!is.na(tgenes)]
tgenes = tgenes[tgenes != '']
tgenes = paste0(tgenes,collapse =',')
inputGOenrich$genes[i] = tgenes
}
writetxt(inputGOenrich,'./out/20230217_aging/Figure2_DEG_GO_tissue/protein/_inputGOenrich_age_DEpro_down.txt',row.names = F,col.names = F)
outids = c()
tpath1 = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_up/Enrichment_heatmap/HeatmapSelectedGO.csv'
tpath2 = './out/20230217_aging/Figure2_DEG_GO_tissue/protein//metascape_DEpro_down/Enrichment_heatmap/HeatmapSelectedGO.csv'
selectGO = readxl::read_xlsx(path = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_up/metascape_result.xlsx',
sheet = 'Enrichment')
go1 = file2frame(tpath1,sep = ',',header = T,row.names =2)
rownames(go1) = paste0(go1$GO,'-',rownames(go1))
go1term = go1$GO
tpath1a = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_up/Enrichment_GO/GO_AllLists.csv'
go1qval = file2frame(tpath1a,sep = ',',header = T)
rownames(go1qval) = paste0(go1qval$GO,'X_LogP_',go1qval$GeneList)
go1 = abs(as.matrix(go1[,-1]))
cname = colnames(go1)
for(i in 1:nrow(go1)){
for(j in 1:ncol(go1)){
tmpname = paste0(go1term[i],cname[j])
go1[i,j] = abs(go1qval[tmpname,]$Log.q.value.)
}
}
go1[is.na(go1)] = 0
go1for_write = cbind(data.frame(GO = rownames(go1)),-go1)
outlist = list()
outlist[[1]] = go1for_write;
outlist[[2]] = selectGO;
tmpnames = sort(unique(go1qval$GeneList))
for(i in 1:length(tmpnames)){
outlist[[i+2]] = go1qval[go1qval$GeneList == tmpnames[[i]],]
}
names(outlist) = c('GOenrichmentALL','selectGO',tmpnames)
openxlsx::write.xlsx(outlist, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S4_GO_tissue_upregulated_protein.xlsx")
selectGO = readxl::read_xlsx(path = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_down/metascape_result.xlsx',sheet = 'Enrichment')
go2 = file2frame(tpath2,sep = ',',header = T,row.names =2)
rownames(go2) = paste0(go2$GO,'-',rownames(go2))
go2term = go2$GO
tpath2a = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_down/Enrichment_GO/GO_AllLists.csv'
go2qval = file2frame(tpath2a,sep = ',',header = T)
rownames(go2qval) = paste0(go2qval$GO,'X_LogP_',go2qval$GeneList)
go2 = abs(as.matrix(go2[,-1]))
cname = colnames(go2)
for(i in 1:nrow(go2)){
for(j in 1:ncol(go2)){
tmpname = paste0(go2term[i],cname[j])
go2[i,j] = -abs(go2qval[tmpname,]$Log.q.value.)
}
}
go2[is.na(go2)] = 0
go1for_write = cbind(data.frame(GO = rownames(go2)),go2)
outlist = list()
outlist[[1]] = go1for_write;
outlist[[2]] = selectGO;
tmpnames = sort(unique(go2qval$GeneList))
for(i in 1:length(tmpnames)){
outlist[[i+2]] = go2qval[go2qval$GeneList == tmpnames[[i]],]
}
names(outlist) = c('GOenrichmentALL','selectGO',tmpnames)
openxlsx::write.xlsx(outlist, file = "./out/20230217_aging/Figure2_DEG_GO_tissue/Data S5_GO_tissue_downregulated_protein.xlsx")
thisgo = rbind(go1,go2)
thisgo.matrix = as.matrix(thisgo)
colnames(thisgo.matrix) = capitalize(gsub('X_LogP_','',colnames(thisgo.matrix)))
thisgo.matrix[abs(thisgo.matrix) < -log10(0.05)] = 0
writetxt(thisgo.matrix,'./out/20230217_aging/Figure2_DEG_GO_tissue/Figure2_heatmap_GOenrichment_metascape_DEpro_eachtissue_fdr.txt',row.names = T)
thisgo.matrix[thisgo.matrix > 4] = 4
thisgo.matrix[thisgo.matrix < -4] = -4
rownames(thisgo.matrix)[rownames(thisgo.matrix) == 'R-HSA-9716542-Signaling by Rho GTPases, Miro GTPases and RHOBTB3'] = 'R-HSA-9716542-Signaling by Rho GTPases'
rownames(thisgo.matrix)[rownames(thisgo.matrix) == 'R-HSA-1428517-The citric acid (TCA) cycle and respiratory electron transport'] = 'R-HSA-1428517-TCA cycle and respiratory electron transport'
graphics.off()
pheatmap::pheatmap(thisgo.matrix,cluster_rows = T,cluster_cols = T,
main ='Enrichment of proteome',
fontsize_row = 11,fontsize_col = 11,fontsize = 14,
treeheight_row = 20,treeheight_col = 20,legend = T,
color=colorRampPalette(c('#3B4992','gray99','#BB0021FF'))(50),
height = 10,width = 10
#file ="./out/20230217_aging/Figure2_DEG_GO_tissue/Figure2_heatmap_GOenrichment_metascape_DEpro_eachtissue.v1_aaas.pdf"
)
7 Figure 3 clustering aging type
7.1 mfuzz promet tissues
tissues = names(pro.tissues)
promet.tissues = list()
promet.tissues.info = list()
promet.tissues.Z = list()
mfuzz.promet.tissues = list()
promet.mstd.eset = list()
for(i in 1:length(tissues)){
#for(i in 1:1){
tt = tissues[i]
thispro = pro.tissues[[tt]]
thispro = delete_dup_genes_forprotein(thispro,pro.tissues.header[[tt]])
thispro.header = pro.tissues.header[[tt]]
thispro = thispro[rowSums(is.na(thispro)) < 1/3*ncol(thispro), ]
thismet = met.tissues[[tt]]
rownames(thismet) = paste0('met_',rownames(thismet))
thispromet = rbind2(thispro,thismet)
thisinfo = pro.tissues.info[[tt]]
thisinfo = thisinfo[colnames(thispromet),]
thispromet.median = t(aggregate(t(thispromet), by=list(thisinfo$stage),
FUN=median, na.rm = T))
mstd = standardise_matrix(thispromet.median)
promet.tissues.Z[[i]] = mstd
promet.tissues[[i]] = thispromet
promet.tissues.info[[i]] = thisinfo
mstd.v = mstd[rowSums(is.na(mstd)) == 0,]
mstd.eset = new("ExpressionSet",exprs = mstd.v)
promet.mstd.eset[[i]] = mstd.eset
mfuzz.promet.tissues[[i]] = mfuzz(mstd.eset, c = 8,m = 1.5)
}
names(mfuzz.promet.tissues) = tissues
names(promet.tissues) = tissues
names(promet.tissues.Z) = tissues
names(promet.tissues.info) = tissues
names(promet.mstd.eset) = tissues
## not run
tissues = names(promet.tissues.mstd.eset)
for(i in 1:length(mfuzz.promet.tissues)){
tpath = paste0('./out/20210428_aging/promet/tissues/',tissues[i])
dir.create(tpath)
pdf(paste0(tpath,'/mfuzz_plot_8A_promet.pdf'),width = 8,height = 4)
mfuzz.plot(promet.tissues.mstd.eset[[i]],mfuzz.promet.tissues[[i]],
mfrow=c(2,4),time.labels = c("Juvenile", "Young_adult","Middle_aged",
"Elderly"),new.window = F)
dev.off()
}
promet.tissues.mstd.eset = promet.mstd.eset
promet.tissues.Z.t = list()
for(i in 1:length(promet.tissues.Z)){
tmp = as.data.frame(t(promet.tissues.Z[[i]]))
promet.tissues.Z.t[[i]] = tmp
}
names(promet.tissues.Z.t) = names(promet.tissues.Z)
promet.whole.Z = t(as.matrix(as.data.frame(rbindlist(promet.tissues.Z.t,fill = T))))
colnames(promet.whole.Z) = paste0(rep(names(promet.tissues.Z),each = 4),'_',rep(1:4,times = 30))
promet.whole.Z.info = data.frame(tissue = rep(names(promet.tissues.Z),each = 4),
stringsAsFactors = F,
stage = rep(1:4,times = 30))
rownames(promet.whole.Z.info) = colnames(promet.whole.Z)
7.2 Figure 3A whole body mfuzz
#this is used to reproduce the figures in manuscript
load('./out/promet_outdata.Rdata')
promet.whole.Z.v = promet.whole.Z[rowSums(is.na(promet.whole.Z)) < 0.5*120,]
promet.whole.Z.mean = t(aggregate(t(promet.whole.Z.v),
by=list(promet.whole.Z.info$stage), FUN=mean, na.rm = T))
promet.whole.Z.mean =promet.whole.Z.mean[-1,]
promet.whole.Z.mean.v = promet.whole.Z.mean[rowSums(is.na(promet.whole.Z.mean)) == 0,]
promet.whole.Z.mean.eset = new("ExpressionSet",exprs = promet.whole.Z.mean.v)
dim(promet.whole.Z.mean.v)
## [1] 5331 4
sum(substr(rownames(promet.whole.Z.mean.v),1,4) == 'met_')
## [1] 1221
#mfuzz.promet.whole = mfuzz(promet.norm.eset.stand, c = 8,m = 1.5)
mfuzz.plot(promet.whole.Z.mean.eset,mfuzz.promet.whole,mfrow=c(2,4),
time.labels = c("Juvenile", "Young_adult","Middle_aged",
"Elderly"),new.window = F)
7.3 Figure 3B Go enrichment
outids = c()
tpath = paste0('./out/20230217_aging/Figure3_trajactory_analysis/cluster1-8_metascape/metascape_result_curated.xlsx')
my_data <- read_excel(tpath, sheet = "Enrichment")
idx = regexpr('Summary',my_data$GroupID) >0
outids = c(outids,my_data$Term[idx])
outids
## [1] "R-HSA-2262752" "GO:0006412" "GO:0045055" "R-HSA-199991"
## [5] "GO:0006163" "GO:0072594" "GO:0006091" "R-HSA-72203"
## [9] "WP3888" "GO:0022411" "GO:0007005" "GO:0006520"
## [13] "GO:0010256" "GO:0097435" "R-HSA-1280215" "GO:0051640"
## [17] "GO:0060627" "GO:1903827" "GO:0006914" "ko04144"
tpath = paste0('./out/20230217_aging/Figure3_trajactory_analysis/cluster1-8_metascape/Enrichment_GO/','GO_membership.csv')
thisgo = file2frame(tpath,sep = ',')
thisgo= thisgo[!duplicated(thisgo$GO),]
xid = is.element(thisgo$GO,outids)
thisgo = thisgo[xid,]
#rownames(thisgo) = paste0(thisgo$GO,':',thisgo$Description)
rownames(thisgo) = capitalize(paste0(thisgo$Description))
thisgo.matrix = as.matrix(thisgo[,substr(colnames(thisgo),1,6)== 'X_LogP'])
writetxt(thisgo.matrix,'./out/20230217_aging/Figure3_trajactory_analysis/Table S1_GOenrichment_metascape_promet_clusters.txt',row.names = T)
thisgo.matrix[thisgo.matrix > -3] = 0
colnames(thisgo.matrix) = capitalize(gsub('X_LogP_','',colnames(thisgo.matrix)))
colnames(thisgo.matrix) = gsub('Cluster','C',colnames(thisgo.matrix))
pheatmap::pheatmap(thisgo.matrix,cluster_rows = T,cluster_cols = T,
main = 'Enrichment analysis of each cluster',
fontsize_row = 9,fontsize_col = 10,fontsize = 10,
treeheight_row = 20,treeheight_col = 20,legend = T,
color=colorRampPalette(c('#3C5488FF','gray95','gray95'))(7),
breaks = c(-32,-24,-16,-8,0),
#file ="./out/20230217_aging/Figure3_trajactory_analysis/Figure3B_GOenrichment_metascape_prometA.pdf",
height = 3.5,width = 5.5
)
7.4 mfuzz for each tissue
7.4.1 get data
# construct data
tissue_trajectory = data.frame()
tissue_names = names(promet.tissues)
for(k in 1:8){
tmpGeneList = names(mfuzz.promet.whole$cluster)[mfuzz.promet.whole$cluster == k]
for(i in 1:length(promet.tissues)){
M1 = promet.tissues[[i]]
metadata.tissue = promet.tissues.info[[i]]
idgene1 = intersect(tmpGeneList,rownames(M1))
cc = repmat(as.matrix(rowMedians(M1[,metadata.tissue$stage == '1'],na.rm = T)),1,ncol(M1))
tsd = repmat(as.matrix(apply(M1,1,sd,na.rm = T)),1,ncol(M1))
M1.Z = (M1 - cc)/tsd
M1.Z.v = M1.Z[idgene1,]
M1.Z.v.mean = t(aggregate(t(M1.Z.v),
by=list(metadata.tissue$stage), FUN=median, na.rm = T))
M1.Z.v.mean = M1.Z.v.mean[-1,]
mean_x = colMeans(M1.Z.v.mean,na.rm =T)
mean_x = mean_x - mean_x[1]
tmpdata2 = data.frame(expr = mean_x,stringsAsFactors = F,
stage = c(1,2,3,4),
group = factor(c('Juvenile','Young','Middle_aged','Elderly'),
level = c('Juvenile','Young','Middle_aged','Elderly')),
cluster = rep(k,4),
tissue = rep(tissue_names[i],4)
)
if (nrow(tissue_trajectory) <1){
tissue_trajectory = tmpdata2
}else{
tissue_trajectory = rbind(tissue_trajectory,tmpdata2)
}
}
}
tissue_trajectory$tissue_system = tissue.systems[tissue_trajectory$tissue]
tissue_trajectory$tissue_system_color = tissue.color[tissue_trajectory$tissue_system]
# to matrix
tissue_trajectory_matrix = matrix(0,nrow(tissue_trajectory)/length(tissue_names),length(tissue_names))
for(i in 1:length(tissue_names)){
tmptr = tissue_trajectory[tissue_trajectory$tissue == tissue_names[i],]
tissue_trajectory_matrix[,i] = tmptr$expr
}
rownames(tissue_trajectory_matrix) = paste(tmptr$group,tmptr$cluster,sep = '_C')
colnames(tissue_trajectory_matrix) = tissue_names
tissue_tr_plot = list()
for (i in 1:8){
idx = tissue_trajectory$cluster == i
tissue_tr_plot [[i]] = ggplot(tissue_trajectory[idx,],
aes(x= group,y = expr,group =tissue)) +
geom_line(aes(color = tissue_system),alpha = 0.3 ,
position = position_dodge(0.2),size = 2)+ scale_color_aaas()+
lghplot.addtheme(hjust = 1,size = 14)+ ggtitle(paste0('Cluster ',i))+
theme(axis.line = element_line(size = 1.2))+ xlab('')+ ylab('')+
scale_y_continuous(labels = scales::comma_format(accuracy =0.1))
}
#pdf(file = "./out/20230217_aging/Figure3_trajactory_analysis/
# Figure3c_trajectory_for_each_tissue_aaas.pdf",height = 9,width = 8)
grid.arrange(arrangeGrob(grobs = tissue_tr_plot,ncol = 3,heights = c(4,4,4),
top = textGrob('Average molecular aging trajectory in 30 tissues',
gp=gpar(fontface="bold", fontsize=20)),
bottom=textGrob('Tissues', gp=gpar(fontface="bold", fontsize=18)),
left = textGrob('Z values', gp=gpar(fontface="bold", fontsize=18),rot=90)))
#dev.off()
8 Figure 4 cluster by MAA
8.1 Figure 4A
promet.norm.eset.stand.matrix = as.matrix(promet.whole.Z.mean.eset)
variability = rep(0,8)
for(j in 1:8){
tgene = names(mfuzz.promet.whole$cluster)[mfuzz.promet.whole$cluster == j]
tgene = intersect(tgene,rownames(promet.norm.eset.stand.matrix))
bx = t(t(promet.norm.eset.stand.matrix[tgene,]) - mfuzz.promet.whole$centers[j,])
variability[j] = mean(sqrt(rowSums(bx^2)/(ncol(bx)-1)))
}
tmpdata = data.frame(amplitude = as.vector(mfuzz.promet.whole$centers[,4]-
mfuzz.promet.whole$centers[,1]),
stringsAsFactors = F,
variability = variability,
class = paste0('Cluster ',1:8))
#pdf(file ='./out/20230217_aging/Figure4_MAA_analysis/
# FigureS4_cluster_amplitude_variability_8_promet.pdf',width = 7,height = 5)
ggplot(tmpdata,aes(variability,amplitude)) + geom_point(size = 4) + lghplot.addtheme(size = 18)+
geom_text_repel(aes(label = class),size = 7)+theme(axis.line = element_line(size = 1.0))+
ggtitle('Figure S4 cluster_amplitude_variability')
#dev.off()
tissues = names(promet.tissues.Z)
clusterdist = matrix(0,length(tissues),8)
clusteramplitude = matrix(0,length(tissues),8)
clusteramplitude_xx = matrix(0,length(tissues),8)
for(i in 1:length(tissues)){
mstd = promet.tissues.Z[[i]]
for(j in 1:8){
tgene = names(mfuzz.promet.whole$cluster)[mfuzz.promet.whole$cluster == j]
tgene = intersect(tgene,rownames(mstd))
bx = t(t(mstd[tgene,]) - mfuzz.promet.whole$centers[j,])
clusterdist[i,j] = mean(sqrt(rowSums(bx^2)/(ncol(bx)-1)),na.rm = T)
if(length(tgene) < 2){
clusteramplitude[i,j] = NA
next;
}
tamp = colMeans(mstd[tgene,],na.rm = T)
clusteramplitude[i,j] = abs(tamp[4]-tamp[1])
clusteramplitude_xx[i,j] = tamp[4]-tamp[1]
}
}
rownames(clusteramplitude) = tissues
rownames(clusteramplitude_xx) = tissues
rownames(clusterdist) = tissues
px = list()
for(j in 1:8){
tmpdata = data.frame(amplitude = clusteramplitude_xx[,j],
stringsAsFactors = F,
variability = clusterdist[,j],
class = tissues,
tissue.systems = tissue.systems)
px[[j]] = ggplot(tmpdata,aes(variability,amplitude,color = tissue.systems)) +
geom_point(size = 6) + lghplot.addtheme(size = 14)+
scale_color_aaas()+scale_fill_aaas()+
geom_text_repel(aes(label = class),size = 3,box.padding = 0.5)+
theme(axis.line = element_line(size = 1.5))+ ggtitle(paste0('Cluster ',j))+
xlab('') + ylab('')
}
#pdf(file = "./out/20230217_aging/Figure4_MAA_analysis/
# FigureXS_cluster8_tissues_amplitude_variability_allV2.pdf",height = 9,width = 16)
grid.arrange(arrangeGrob(grobs = px,ncol = 4,margin = c(1,1,1,1),
top = textGrob('Molecular alteration amplitude of each cluster in 30 tissues',
gp=gpar(fontface="bold", fontsize=24)),
bottom=textGrob('Molecular alteration variability',
gp=gpar(fontface="bold", fontsize=24)),
#bottom = 'mRNA expression(log2 FPKM)',
left = textGrob('Molecular alteration amplitude(MAA)',
gp=gpar(fontface="bold", fontsize=24),rot=90)))
#dev.off()
8.3 Figure 4C
dendrow = as.dendrogram(dent.MAA$tree_row)
labelColors = c('#3C5488FF','#E64B35FF')
clusMember = cutree(dent.MAA$tree_row,2)
# function to get color labels
colLab <- function(n) {
if (is.leaf(n)) {
a <- attributes(n)
labCol <- labelColors[clusMember[which(names(clusMember) == a$label)]]
attr(n, "nodePar") <- c(a$nodePar, lab.col = labCol)
}
n
}
clusDendro = dendrapply(dendrow, colLab)
#pdf('./out/20230217_aging/Figure4_MAA_analysis/
# FigureX_dendrogam_by_MAA.pdf',width = 7,height = 5)
plot(clusDendro, main = "Molecular trajactory dendgrogram by MAA")
#dev.off()
mean_clusteramplitude_xx = rowMeans(clusteramplitude_xx)
names(mean_clusteramplitude_xx) = rownames(clusteramplitude_xx)
tissues_pro = rownames(clusteramplitude_xx)
dendrow = as.dendrogram(dent.MAA$tree_row)
dendrowsplit = cut(dendrow,2.5)
class1 = unlist(cut(dendrow,2.5)$lower[[1]])
class2 = unlist(cut(dendrow,2.5)$lower[[2]])
tissueClass.MAA = data.frame(tissue = c(tissues_pro[class1],
tissues_pro[class2]),stringsAsFactors = F,
class = c(rep('Type II',length(class1)),
rep('Type I',length(class2))))
rownames(tissueClass.MAA) = tissueClass.MAA$tissue
tissueClass.MAA$mean_aging_amplitude=mean_clusteramplitude_xx[rownames(tissueClass.MAA)]
tissueClass$class.MAA = tissueClass.MAA[rownames(tissueClass),]$class
tissueClass$mean_aging_amplitude = tissueClass.MAA[rownames(tissueClass),]$mean_aging_amplitude
tissueClass$type = rep('Undefine',nrow(tissueClass))
tissueClass$type[tissueClass$class == 'Type I' & tissueClass$class.MAA == 'Type I'] = 'Type I'
tissueClass$type[tissueClass$class == 'Type II' & tissueClass$class.MAA == 'Type II'] = 'Type II'
tissueClass
## tissue type class class.MAA
## Pancreas Pancreas Type II Type II Type II
## Supramarginal_gyrus Supramarginal_gyrus Type II Type II Type II
## Frontal_pole Frontal_pole Type II Type II Type II
## Liver Liver Type II Type II Type II
## Hippocampus Hippocampus Type II Type II Type II
## Skin_of_back Skin_of_back Type II Type II Type II
## Cecum Cecum Type II Type II Type II
## Lung Lung Type II Type II Type II
## Superior_temporal_gyrus Superior_temporal_gyrus Type II Type II Type II
## Muscle Muscle Type II Type II Type II
## Adrenal_gland Adrenal_gland Type II Type II Type II
## Hypothalamus Hypothalamus Type II Type II Type II
## Heart Heart Type II Type II Type II
## Arteria_cruralis Arteria_cruralis Undefine Type II Type I
## Arteria_carotis Arteria_carotis Undefine Type II Type I
## Adipose Adipose Type I Type I Type I
## Thyroid_gland Thyroid_gland Type I Type I Type I
## Facial_skin Facial_skin Undefine Type I Type II
## Femoral_vein Femoral_vein Undefine Type I Type II
## Ileocecum Ileocecum Type I Type I Type I
## Thymus Thymus Type I Type I Type I
## Aortic_arch Aortic_arch Type I Type I Type I
## Spleen Spleen Type I Type I Type I
## Kidney Kidney Type I Type I Type I
## Ovary Ovary Type I Type I Type I
## Fallopian_tube Fallopian_tube Type I Type I Type I
## Stomach Stomach Type I Type I Type I
## Duodenum Duodenum Type I Type I Type I
## Pituitary Pituitary Type I Type I Type I
## Uterus Uterus Type I Type I Type I
## mean_aging_amplitude
## Pancreas -0.075748775
## Supramarginal_gyrus 0.231195111
## Frontal_pole 0.186849533
## Liver 0.210951149
## Hippocampus 0.044697652
## Skin_of_back 0.153227025
## Cecum 0.117593492
## Lung 0.062388571
## Superior_temporal_gyrus 0.151674043
## Muscle 0.127531951
## Adrenal_gland 0.390895560
## Hypothalamus 0.312324819
## Heart 0.204160794
## Arteria_cruralis -0.302318429
## Arteria_carotis -0.015271359
## Adipose -0.650050564
## Thyroid_gland -0.233352138
## Facial_skin -0.068819684
## Femoral_vein 0.002050297
## Ileocecum -0.157174219
## Thymus -0.293446352
## Aortic_arch -0.341411484
## Spleen -0.253614804
## Kidney -0.119119447
## Ovary -0.148487802
## Fallopian_tube -0.242323566
## Stomach -0.151051049
## Duodenum -0.334290753
## Pituitary -0.074656135
## Uterus -0.150139727
require(scatterplot3d)
#typecolor = '#3B4992','gray99','#EE0000FF',
#pdf('./out/20230217_aging/Figure4_MAA_analysis/Figure4C_3d_scatterplotAA.pdf')
typecolor = tissueClass[rownames(clusteramplitude_xx),]$color
s3d <- scatterplot3d(clusteramplitude_xx[,8], clusteramplitude_xx[,2], clusteramplitude_xx[,4], grid = T,size = 1, cex.lab = 2, cex.symbols = 4, color=typecolor,
pch=19,
scale.y=.75,
main="Aging amplitude",
xlab='C8',
ylab='C2 ',
zlab='C4')
s3d1 <- scatterplot3d(clusteramplitude_xx[,8], clusteramplitude_xx[,2], clusteramplitude_xx[,4], grid = T,size = 1, cex.lab = 2, cex.symbols = 4, # x y and z axis
color=typecolor,
pch=19,
scale.y=0.75,
scale.Z=1.25,
main="Aging amplitude",
xlab='C8',
ylab='C2 ',
zlab='C4',
text(s3d$xyz.convert(clusteramplitude_xx[,c(8,2,4)]), labels = rownames(clusteramplitude_xx),
cex= 1))
#dev.off()
9 Figure 5 experimental validation
# Figure 5A all changed protein
tx = colSums(abs(Aging_pro_sigall_matrix),na.rm = T)
yy1 = tissueClass[names(tx),]
#pdf('./out/20230217_aging/Figure5_markers_HE/Figure_S5_allchanged_proteins.pdf')
ggplot(,aes(x = yy1$type,y = tx,color = yy1$type,fill = yy1$type)) +#geom_violin()+
geom_boxplot(outlier.size = -1,alpha = 0.1,size = 0.8)+
geom_jitter(size = 5,width = 0.3)+ theme_classic() +
ggtitle('# of changed proteins')+
lghplot.addtheme()+ scale_color_manual(values=c('#EE0000','#3B4992','gray10'))+
xlab('Tissue Type')+ylab('Number of changed proteins')
#dev.off()
wilcox.test(tx~yy1$type,subset= yy1$type != 'Undefine')
##
## Wilcoxon rank sum exact test
##
## data: tx by yy1$type
## W = 143, p-value = 0.001914
## alternative hypothesis: true location shift is not equal to 0
tx = colSums(abs(Aging_pro_sigall_matrix),na.rm = T)/
colSums(!is.na(Aging_pro_sigall_matrix))
yy1 = tissueClass[names(tx),]
#pdf('./out/20230217_aging/Figure5_markers_HE/Figure4B_allchanged_proteins_percentageV2.pdf')
ggplot(,aes(x = yy1$type,y = tx*100,color = yy1$type,fill = yy1$type)) +#geom_violin()+
geom_boxplot(outlier.size = -1,alpha = 0.1,size = 0.8)+
geom_jitter(size = 5,width = 0.3)+ theme_classic() +
ggtitle('% of changed proteins')+
lghplot.addtheme()+ scale_color_manual(values=c('#EE0000','#3B4992','gray10'))+
xlab('Tissue Type')+ylab('Percentage of changed proteins(%)')
#dev.off()
wilcox.test(tx~yy1$type,subset= yy1$type != 'Undefine')
##
## Wilcoxon rank sum exact test
##
## data: tx by yy1$type
## W = 136, p-value = 0.007244
## alternative hypothesis: true location shift is not equal to 0
9.1 Figure 5B
p16 = file2frame('./data/P16_tissues_v20230419.txt')
utissue = unique(p16$Tissue)
tpvalues = data.frame(tissue = utissue,
p_Juvenile_vs_Elderly = rep(1,length(utissue)),
p_Young_vs_Elderly = rep(1,length(utissue)))
breaks_A = c(0.01,0.08,0.08,0.08,0.08,0.01,0.01,0.01)
breaks_B = c(3,3,3,3,2,2,0.06,3)
for(i in 1:length(utissue)){
tmpdata = p16[p16$Tissue == utissue[i],]
tmpdata$Area.Density.log2 = tmpdata$Area.Density*1e3#log2(tmpdata$Area.Density*1e3)
tmpdata$class = factor(x = tmpdata$class,levels = c('Juvenile','Young','Elderly'))
tmpp1 = t.test(Area.Density~class,data = tmpdata,subset = tmpdata$class != 'Young')$p.value
tmpp2 = t.test(Area.Density~class,data = tmpdata,subset = tmpdata$class != 'Juvenile')$p.value
tpvalues$p_Juvenile_vs_Elderly[i] = tmpp1
tpvalues$p_Young_vs_Elderly[i] = tmpp2
df2 = data.frame(tmean = aggregate(tmpdata$Area.Density.log2,
by=list(tmpdata$class), FUN=mean, na.rm = T)$x,
tsd = aggregate(tmpdata$Area.Density.log2, by=list(tmpdata$class),
FUN=sd, na.rm = T)$x,
class = factor(x= c('Juvenile','Young','Elderly'),
levels = c('Juvenile','Young','Elderly'))
)
thepath = paste0('./out/20230217_aging/Figure5_markers_HE/Figure5x_p16_stats', utissue[i],'.pdf')
pdf(thepath)
p = ggplot(df2, aes(x = class, y = tmean, fill = class)) +
geom_bar(stat = 'identity',color = 'black',position = position_dodge(),
alpha = 0.5,size = 1.5)+
geom_errorbar(aes(ymin = tmean,ymax = tmean+tsd),width = .3,size = 1.5)+
lghplot.addthemeA(size = 28,sizex = 26,sizey = 26)+
scale_fill_manual(values = c('#008B45FF','#3B4992FF','#EE0000FF'))+
theme(axis.line = element_line(size = 1.2))+ xlab('')+
ylab(bquote('Area Density ' ~ italic(x10) ^italic(3)))+
#ylab('Log10 Area Density x 1e-7')+
ggtitle(utissue[i])
#if(i <=8){
# p = p+ scale_y_break(c(breaks_A[i],breaks_B[i]), scales = "free",space = 0.2)
#}
print(p)
dev.off()
}
tpvalues
## tissue p_Juvenile_vs_Elderly p_Young_vs_Elderly
## 1 Aortic_arch 4.546295e-03 4.553667e-03
## 2 Spleen 2.928558e-03 2.941005e-03
## 3 Kidney 3.607333e-04 3.618602e-04
## 4 Ovary 2.399612e-03 2.487708e-03
## 5 Thymus 1.122249e-02 1.151782e-02
## 6 Uterus 1.052903e-04 1.057245e-04
## 7 Thyroid_gland 6.633705e-04 6.751286e-04
## 8 Stomach 3.601685e-05 3.616963e-05
## 9 Pancreas 2.616949e-02 5.453556e-03
## 10 Skin_of_back 2.446174e-02 4.146347e-02
## 11 Lung 4.895734e-03 6.619731e-03
## 12 Liver 2.051557e-01 2.593903e-01
## 13 Muscle 3.401860e-01 6.380951e-01
9.2 Figure 5C
p21 = file2frame('./data/P21_tissues_v20230419.txt')
utissue = unique(p21$Tissue)
tpvalues = data.frame(tissue = utissue,
p_Juvenile_vs_Elderly = rep(1,length(utissue)),
p_Young_vs_Elderly = rep(1,length(utissue)))
breaks_A = c(0.05,0.4,0.4,0.4,0.1,0.4,0.4,0.05)
breaks_B = c(1,4,4,4,2,4,4,1)
for(i in 1:length(utissue)){
tmpdata = p21[p21$Tissue == utissue[i],]
tmpdata$Area.Density.log2 = tmpdata$Area.Density*1e3#log2(tmpdata$Area.Density*1e3)
tmpdata$class = factor(x = tmpdata$class,levels = c('Juvenile','Young','Elderly'))
tmpp1 = t.test(Area.Density~class,data = tmpdata,subset = tmpdata$class != 'Young')$p.value
tmpp2 = t.test(Area.Density~class,data = tmpdata,subset = tmpdata$class != 'Juvenile')$p.value
tpvalues$p_Juvenile_vs_Elderly[i] = tmpp1
tpvalues$p_Young_vs_Elderly[i] = tmpp2
df2 = data.frame(tmean = aggregate(tmpdata$Area.Density.log2,
by=list(tmpdata$class), FUN=mean, na.rm = T)$x,
tsd = aggregate(tmpdata$Area.Density.log2, by=list(tmpdata$class),
FUN=sd, na.rm = T)$x,
class = factor(x= c('Juvenile','Young','Elderly'),
levels = c('Juvenile','Young','Elderly'))
)
thepath = paste0('./out/20230217_aging/Figure5_markers_HE/Figure4x_p21_stats', utissue[i],'.pdf')
pdf(thepath)
p = ggplot(df2, aes(x = class, y = tmean, fill = class)) +
geom_bar(stat = 'identity',color = 'black',position = position_dodge(),
alpha = 0.5,size = 1.5)+
geom_errorbar(aes(ymin = tmean,ymax = tmean+tsd),width = .3,size = 1.5)+
lghplot.addthemeA(size = 28,sizex = 26,sizey = 26)+
scale_fill_manual(values = c('#008B45FF','#3B4992FF','#EE0000FF'))+
theme(axis.line = element_line(size = 1.2))+ xlab('')+
ylab(bquote('Area Density ' ~ italic(x10) ^italic(3)))+
#ylab('Log10 Area Density x 1e-7')+
ggtitle(utissue[i])
#if(i <=8){
# p = p+ scale_y_break(c(breaks_A[i],breaks_B[i]), scales = "free",space = 0.2)
#}
print(p)
dev.off()
}
tpvalues
## tissue p_Juvenile_vs_Elderly p_Young_vs_Elderly
## 1 Aortic_arch 3.371849e-03 3.490217e-03
## 2 Spleen 6.244567e-04 6.370290e-04
## 3 Kidney 1.052679e-02 1.070538e-02
## 4 Ovary 3.484938e-04 3.559713e-04
## 5 Thymus 4.256548e-02 4.462434e-02
## 6 Uterus 2.571205e-02 2.611076e-02
## 7 Thyroid_gland 1.553183e-02 1.686152e-02
## 8 Stomach 4.794586e-05 4.974984e-05
## 9 Pancreas 2.700608e-01 1.668549e-01
## 10 Skin_of_back 5.553932e-02 2.923295e-01
## 11 Lung 6.940087e-04 6.950478e-04
## 12 Liver 9.815520e-01 9.972549e-01
## 13 Muscle 3.620859e-01 7.890578e-01
9.3 Figure 5D
Cellcounts = file2frame('./data/cell_counts_v20230407.txt')
utissue = unique(Cellcounts$Tissue)
tpvalues = data.frame(tissue = utissue,
p_Juvenile_vs_Elderly = rep(1,length(utissue)),
p_Young_vs_Elderly = rep(1,length(utissue)))
for(i in 1:length(utissue)){
tmpdata = Cellcounts[Cellcounts$Tissue == utissue[i],]
tmpdata$class = factor(x = tmpdata$class,
levels = c('Juvenile','Young','Elderly'))
tmpp1 = t.test(number~class,data = tmpdata,subset = tmpdata$class != 'Young')$p.value
tmpp2 = t.test(number~class,data = tmpdata,subset = tmpdata$class != 'Juvenile')$p.value
tpvalues$p_Juvenile_vs_Elderly[i] = tmpp1
tpvalues$p_Young_vs_Elderly[i] = tmpp2
df2 = data.frame(tmean = aggregate(tmpdata$number, by=list(tmpdata$class),
FUN=mean, na.rm = T)$x,
tsd = aggregate(tmpdata$number, by=list(tmpdata$class), FUN=sd, na.rm = T)$x,
class = factor(x= c('Juvenile','Young','Elderly'),
levels = c('Juvenile','Young','Elderly'))
)
thepath = paste0('./out/20230217_aging/Figure5_markers_HE/Figure4x_HE_cellcounts_', utissue[i],'_v1.pdf')
pdf(thepath)
p = ggplot(df2, aes(x = class, y = tmean, fill = class)) +
geom_bar(stat = 'identity',color = 'black',position = position_dodge(),alpha = 0.5,size = 1.5)+
geom_errorbar(aes(ymin = tmean,ymax = tmean+tsd),width = .3,size = 1.5)+
lghplot.addthemeA(size = 28,sizex = 26,sizey = 26)+
scale_fill_manual(values = c('#008B45FF','#3B4992FF','#EE0000FF'))+
theme(axis.line = element_line(size = 1.2))+ xlab('')+
ylab('Number of Parenchymal Cells')+ggtitle(utissue[i])
print(p)
dev.off()
}
tpvalues
## tissue p_Juvenile_vs_Elderly p_Young_vs_Elderly
## 1 Thymus 3.842297e-07 2.479669e-04
## 2 Stomach 4.735346e-08 4.458945e-06
## 3 Aortic_arch 8.038428e-07 2.041040e-07
## 4 Ovary 6.576220e-09 6.727973e-05
## 5 Spleen 1.724836e-06 4.997698e-01
## 6 Thyroid_gland 3.814005e-06 9.571183e-02
## 7 Kidney 2.734793e-01 3.090672e-01
10 Figure 6: Translation efficiency
10.1 ratio data construction
tissues = names(promet.tissues.Z)
clusterdist = matrix(0,length(tissues),8)
clusteramplitude = matrix(0,length(tissues),8)
clusteramplitude_xx = matrix(0,length(tissues),8)
for(i in 1:length(tissues)){
mstd = promet.tissues.Z[[i]]
for(j in 1:8){
tgene = names(mfuzz.promet.whole$cluster)[mfuzz.promet.whole$cluster == j]
tgene = intersect(tgene,rownames(mstd))
bx = t(t(mstd[tgene,]) - mfuzz.promet.whole$centers[j,])
clusterdist[i,j] = mean(sqrt(rowSums(bx^2)/(ncol(bx)-1)),na.rm = T)
if(length(tgene) < 2){
clusteramplitude[i,j] = NA
next;
}
tamp = colMeans(mstd[tgene,],na.rm = T)
clusteramplitude[i,j] = abs(tamp[4]-tamp[1])
clusteramplitude_xx[i,j] = tamp[4]-tamp[1]
}
}
rownames(clusteramplitude) = tissues
rownames(clusteramplitude_xx) = tissues
rownames(clusterdist) = tissues
ratio.tissues = list()
ratio.tissues.info = list()
pro.tissues.forRatio = list()
mrna.tissues.forRatio = list()
tissues = names(pro.tissues)
overlaptissues = intersect(names(pro.tissues),names(mrna.tissues))
for( i in 1:length(overlaptissues)){
this_tissue = overlaptissues[i]
thispro = delete_dup_genes_forprotein(pro.tissues[[this_tissue]],
pro.tissues.header[[this_tissue]])
thismrna = mrna.tissues[[this_tissue]]
vcol = intersect(colnames(thispro),colnames(thismrna))
vrow = intersect(rownames(thispro),rownames(thismrna))
thispro = thispro[vrow,vcol]
thismrna = thismrna[vrow,vcol]
pro.tissues.forRatio[[i]] = thispro
mrna.tissues.forRatio[[i]] = thismrna
ratio.tissues[[i]] = thispro - thismrna
ratio.tissues.info[[i]] = pro.tissues.info[[i]][vcol,]
}
names(ratio.tissues) = overlaptissues
names(pro.tissues.forRatio) = overlaptissues
names(mrna.tissues.forRatio) = overlaptissues
names(ratio.tissues.info) = overlaptissues
ratio_amplitude = list()
pro_amplitude = list()
mrna_amplitude = list()
n = length(ratio.tissues)
ratio_out = data.frame(meanChangeRatio = rep(0,n),stringsAsFactors = F,
fc_up_down = rep(0,n),
tissues = names(ratio.tissues),
meanChangeproZ = rep(0,n),
fc_up_down_pro = rep(0,n),
meanChangemrnaZ = rep(0,n),
fc_up_down_mrna = rep(0,n))
for (i in 1:length(ratio.tissues)){
bx = ratio.tissues[[i]]
bx.info = ratio.tissues.info[[i]]
id = bx.info$stage < 5
bx = bx[,id]
bx.info = bx.info[id,]
cx = t(aggregate(t(bx), by=list(bx.info$stage), FUN=mean, na.rm = T))
cx = as.matrix(standardise_1(new("ExpressionSet",exprs = cx)))
cx = cx[-1,]
ee = cx[,4]-cx[,1]
ratio_amplitude[[i]] = ee;
#hist(ee,30)
ratio_out$meanChangeRatio[i] = mean(ee,na.rm = T)
ratio_out$fc_up_down[i] = log2(sum(ee > 0,na.rm = T)/sum(ee < -0,na.rm = T))
bx = pro.tissues.forRatio[[i]]
bx.info = ratio.tissues.info[[i]]
id = bx.info$stage < 5
bx = bx[,id]
bx.info = bx.info[id,]
cx = t(aggregate(t(bx), by=list(bx.info$stage), FUN=mean, na.rm = T))
cx = as.matrix(standardise_1(new("ExpressionSet",exprs = cx)))
cx = cx[-1,]
ee = cx[,4]-cx[,1]
pro_amplitude[[i]] = ee
#hist(ee,30)
ratio_out$meanChangeproZ[i] = mean(ee,na.rm = T)
ratio_out$fc_up_down_pro[i] = log2(sum(ee > 0,na.rm = T)/sum(ee < -0,na.rm = T))
bx = mrna.tissues.forRatio[[i]]
bx.info = ratio.tissues.info[[i]]
id = bx.info$stage < 5
bx = bx[,id]
bx.info = bx.info[id,]
cx = t(aggregate(t(bx), by=list(bx.info$stage), FUN=mean, na.rm = T))
cx = as.matrix(standardise_1(new("ExpressionSet",exprs = cx)))
cx = cx[-1,]
ee = cx[,4]-cx[,1]
mrna_amplitude[[i]] = ee
#hist(ee,30)
ratio_out$meanChangemrnaZ[i] = mean(ee,na.rm = T)
ratio_out$fc_up_down_mrna[i] = log2(sum(ee > 0,na.rm = T)/sum(ee < -0,na.rm = T))
}
rownames(ratio_out) = names(ratio.tissues)
names(ratio_amplitude) = names(ratio.tissues)
names(pro_amplitude) = names(ratio.tissues)
names(mrna_amplitude) = names(ratio.tissues)
10.2 Figure 6A design
11 Figure 7
11.1 GO data construct
#outids = c('R-HSA-72766','GO:0006413','hsa03010','ko04142','hsa04120','ko03050')
outids = c('R-HSA-72766','hsa03010','hsa04142','hsa04120','hsa03050')
#tpath = paste0('./out/20210428_aging/promet/tissues/metascape/Aging_up_genes_metascape/Enrichment_GO/','GO_membership.csv')
tpath = paste0('./out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_up/Enrichment_GO/','GO_membership.csv')
thisgo = file2frame(tpath,sep = ',')
thisgo= thisgo[!duplicated(thisgo$GO),]
xid = is.element(thisgo$GO,outids)
thisgo = thisgo[xid,]
#rownames(thisgo) = paste0(thisgo$GO,':',thisgo$Description)
rownames(thisgo) = paste0(thisgo$GO)
thisgo.matrix = -as.matrix(thisgo[,substr(colnames(thisgo),1,6)== 'X_LogP'])
# replace qval
tpath1a = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_up/Enrichment_GO/GO_AllLists.csv'
go1qval = file2frame(tpath1a,sep = ',',header = T)
rownames(go1qval) = paste0(go1qval$GO,'X_LogP_',go1qval$GeneList)
go1term = rownames(thisgo.matrix)
cname = colnames(thisgo.matrix)
for(i in 1:nrow(thisgo.matrix)){
for(j in 1:ncol(thisgo.matrix)){
tmpname = paste0(go1term[i],cname[j])
thisgo.matrix[i,j] = abs(go1qval[tmpname,]$Log.q.value.)
}
}
thisgo.matrix[is.na(thisgo.matrix)] = 0
thisgo.matrix[abs(thisgo.matrix) < 3] = 0
colnames(thisgo.matrix) = capitalize(gsub('X_LogP_','',colnames(thisgo.matrix)))
colnames(thisgo.matrix) = gsub('Cluster','C',colnames(thisgo.matrix))
thisgo.matrix.up = thisgo.matrix
tmpname = rownames(thisgo.matrix.up)
#add hsa03050
thisgo.matrix.up = rbind(thisgo.matrix.up,rep(0,30))
rownames(thisgo.matrix.up) = c(tmpname,'hsa03050')
thisgo.matrix.up = thisgo.matrix.up[outids,]
thisgo.matrix.up
## Adipose Adrenal_gland Aortic_arch Arteria_carotis Arteria_cruralis
## R-HSA-72766 0 0 0 0 0
## hsa03010 0 0 0 0 0
## hsa04142 0 0 0 0 0
## hsa04120 0 0 0 0 0
## hsa03050 0 0 0 0 0
## Cecum Duodenum Facial_skin Fallopian_tube Femoral_vein Frontal_pole
## R-HSA-72766 0 0 0 0 4.3 0
## hsa03010 0 0 0 0 0.0 0
## hsa04142 0 0 0 0 0.0 0
## hsa04120 0 0 0 0 0.0 0
## hsa03050 0 0 0 0 0.0 0
## Heart Hippocampus Hypothalamus Ileocecum Kidney Liver Lung Muscle
## R-HSA-72766 0.0 0 0 0 0 0.0 0 0
## hsa03010 0.0 0 0 0 0 0.0 0 0
## hsa04142 3.5 0 3 0 0 3.2 0 0
## hsa04120 0.0 0 0 0 0 0.0 0 0
## hsa03050 0.0 0 0 0 0 0.0 0 0
## Ovary Pancreas Pituitary Skin_of_back Spleen Stomach
## R-HSA-72766 0 0 0 0 0 0
## hsa03010 0 0 0 0 0 0
## hsa04142 0 0 0 0 0 0
## hsa04120 0 0 0 0 0 0
## hsa03050 0 0 0 0 0 0
## Superior_temporal_gyrus Supramarginal_gyrus Thymus Thyroid_gland
## R-HSA-72766 0 0 0 0
## hsa03010 0 0 0 0
## hsa04142 0 0 0 0
## hsa04120 0 0 0 0
## hsa03050 0 0 0 0
## Uterus
## R-HSA-72766 0
## hsa03010 0
## hsa04142 0
## hsa04120 0
## hsa03050 0
outids = c('R-HSA-72766','hsa03010','hsa04142','hsa04120','hsa03050')
outids_des = c('Translation','Ribosome','Lysosome',
'Ubiquitin mediated proteolysis', 'Proteasome')
tpath = paste0('./out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_down/Enrichment_GO/','GO_membership.csv')
thisgo = file2frame(tpath,sep = ',')
thisgo= thisgo[!duplicated(thisgo$GO),]
xid = is.element(thisgo$GO,outids)
thisgo = thisgo[xid,]
rownames(thisgo) = paste0(thisgo$GO)
# replace qval
tpath1a = './out/20230217_aging/Figure2_DEG_GO_tissue/protein/metascape_DEpro_down/Enrichment_GO/GO_AllLists.csv'
go1qval = file2frame(tpath1a,sep = ',',header = T)
rownames(go1qval) = paste0(go1qval$GO,'X_LogP_',go1qval$GeneList)
go1term = rownames(thisgo.matrix)
cname = colnames(thisgo.matrix)
for(i in 1:nrow(thisgo.matrix)){
for(j in 1:ncol(thisgo.matrix)){
tmpname = paste0(go1term[i],cname[j])
thisgo.matrix[i,j] = abs(go1qval[tmpname,]$Log.q.value.)
}
}
thisgo.matrix[is.na(thisgo.matrix)] = 0
thisgo.matrix = as.matrix(thisgo[,substr(colnames(thisgo),1,6)== 'X_LogP'])
thisgo.matrix[abs(thisgo.matrix) < 3 ] = 0
colnames(thisgo.matrix) = capitalize(gsub('X_LogP_','',colnames(thisgo.matrix)))
colnames(thisgo.matrix) = gsub('Cluster','C',colnames(thisgo.matrix))
rownames(thisgo.matrix)
## [1] "R-HSA-72766" "hsa03010" "hsa03050" "hsa04142" "hsa04120"
thisgo.matrix.down = thisgo.matrix
thisgo.matrix.down = thisgo.matrix.down[outids,]
thisgo.matrix.down
## Adipose Adrenal_gland Aortic_arch Arteria_carotis Arteria_cruralis
## R-HSA-72766 -5.2 0 -5.9 0 -4.3
## hsa03010 0.0 0 0.0 0 0.0
## hsa04142 0.0 0 0.0 0 0.0
## hsa04120 0.0 0 0.0 0 0.0
## hsa03050 0.0 0 0.0 0 0.0
## Cecum Duodenum Facial_skin Fallopian_tube Femoral_vein Frontal_pole
## R-HSA-72766 0 0.0 0 -6.9 0 0
## hsa03010 0 0.0 0 -5.0 0 0
## hsa04142 0 0.0 0 0.0 0 0
## hsa04120 0 0.0 0 0.0 0 0
## hsa03050 0 -4.5 0 0.0 0 0
## Heart Hippocampus Hypothalamus Ileocecum Kidney Liver Lung Muscle
## R-HSA-72766 0 0 0 -7.2 0 -6 0 -4.6
## hsa03010 0 0 0 -3.9 0 0 0 0.0
## hsa04142 0 0 0 -3.6 0 0 0 0.0
## hsa04120 0 0 0 0.0 0 0 0 0.0
## hsa03050 0 0 0 0.0 0 0 0 0.0
## Ovary Pancreas Pituitary Skin_of_back Spleen Stomach
## R-HSA-72766 -7.7 0 -11 -3.8 0 -12.0
## hsa03010 0.0 0 -10 0.0 0 -4.7
## hsa04142 0.0 0 0 0.0 0 -5.7
## hsa04120 0.0 0 0 0.0 0 0.0
## hsa03050 0.0 0 0 -3.1 0 -6.7
## Superior_temporal_gyrus Supramarginal_gyrus Thymus Thyroid_gland
## R-HSA-72766 0 0 -34.0 0
## hsa03010 0 0 -7.3 0
## hsa04142 0 0 0.0 0
## hsa04120 0 0 -7.8 0
## hsa03050 0 0 -25.0 0
## Uterus
## R-HSA-72766 -5
## hsa03010 0
## hsa04142 0
## hsa04120 0
## hsa03050 0
11.2 Figure 7A
outids = c('R-HSA-72766','hsa03010','hsa04142','hsa04120','hsa03050')
outids_des = c('Translation','Ribosome','Lysosome',
'Ubiquitin mediated proteolysis', 'Proteasome')
thisgo.matrix.all = thisgo.matrix.up+ thisgo.matrix.down
tmpnames = rownames(tissueClass)
tmpnames = c(tmpnames[tissueClass$type == 'Type II'],tmpnames[tissueClass$type == 'Undefine'],
tmpnames[tissueClass$type == 'Type I'])
thisgo.matrix.all = thisgo.matrix.all[,tmpnames]
thisgo.matrix.all[thisgo.matrix.all >= 3] = 1
thisgo.matrix.all[thisgo.matrix.all <= -3] = -1
tclass = data.frame(class = tissueClass[tmpnames,]$type,row.names = tmpnames)
colnames(tclass) <- c("Tissue_type")
ann_colors = list(
Tissue_type= c('#FD60A7','gray85','#4fffF7')
)
names(ann_colors$Tissue_type) = c('Type I','Undefine', 'Type II')
#sort tissues
tmpname = colnames(thisgo.matrix.all)
rownames(thisgo.matrix.all) = outids_des
pheatmap::pheatmap(thisgo.matrix.all,cluster_rows = F,annotation_colors = ann_colors,
annotation_col = tclass,
main = 'Enrichment in mRNA translation and protein degradation related pathways',
cluster_cols = F,fontsize_row = 9,fontsize_col = 10,
fontsize = 10,treeheight_row = 20,treeheight_col = 20,legend = T,
color=colorRampPalette(c('#3B4992','gray99','#EE0000FF'))(6),
#file ="./out/20230217_aging/Figure6_heatmap_TED_machanism/Figure 6A_enrichment_metascape_promet_all_tissue_TEDs_V6.pdf",
height = 3,width = 10
)
idx = tclass$Tissue_type != 'Undefine'
tmpaa = thisgo.matrix.all[,idx]
tmpclass = tclass$Tissue_type[idx]
rownames(tmpaa)
## [1] "Translation" "Ribosome"
## [3] "Lysosome" "Ubiquitin mediated proteolysis"
## [5] "Proteasome"
cor.test(abs(tmpaa[1,]),(tmpclass == 'Type I')+0,method = 'spearman')
##
## Spearman's rank correlation rho
##
## data: abs(tmpaa[1, ]) and (tmpclass == "Type I") + 0
## S = 1571, p-value = 0.01725
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.46291
cor.test(abs(tmpaa[2,]),(tmpclass == 'Type I')+0,method = 'spearman')
##
## Spearman's rank correlation rho
##
## data: abs(tmpaa[2, ]) and (tmpclass == "Type I") + 0
## S = 1497.7, p-value = 0.01144
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.48795
cor.test(abs(tmpaa[3,]),(tmpclass == 'Type I')+0,method = 'spearman')
##
## Spearman's rank correlation rho
##
## data: abs(tmpaa[3, ]) and (tmpclass == "Type I") + 0
## S = 3210.5, p-value = 0.6353
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.09759001
cor.test(abs(tmpaa[4,]),(tmpclass == 'Type I')+0,method = 'spearman')
##
## Spearman's rank correlation rho
##
## data: abs(tmpaa[4, ]) and (tmpclass == "Type I") + 0
## S = 2340, p-value = 0.3273
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.2
cor.test(abs(tmpaa[5,]),(tmpclass == 'Type I')+0,method = 'spearman')
##
## Spearman's rank correlation rho
##
## data: abs(tmpaa[5, ]) and (tmpclass == "Type I") + 0
## S = 2301.4, p-value = 0.2957
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.2132007
11.3 Figure 7B
require(data.table)
amplitude_matrix = list()
for(i in 1:length(pro_amplitude)){
tmp = matrix(pro_amplitude[[i]],1,length(pro_amplitude[[i]]))
tmp = as.data.frame(tmp)
colnames(tmp) = names(pro_amplitude[[i]])
amplitude_matrix[[i]] = tmp
}
#amplitude_matrix = rbindlist(amplitude_matrix,fill = T)
amplitude_matrix = t(as.matrix(rbindlist(amplitude_matrix,fill = T)))
colnames(amplitude_matrix) = names(pro_amplitude)
vid = rowSums(is.na(amplitude_matrix)) < ncol(amplitude_matrix)/2
amplitude_matrix.v = amplitude_matrix[vid,]
amplitude_matrix.v[is.na(amplitude_matrix.v)] = NA
dim(amplitude_matrix.v)
## [1] 3902 30
tcor = cor(t(amplitude_matrix.v),ratio_out$meanChangeRatio,use = "pairwise")
names(tcor) = rownames(amplitude_matrix.v)
gsea_input = sort(tcor,decreasing = T)
calcp <- function(x,y){
return(cor.test(x,y,use = "pairwise")$p.value)
}
tpvalue = apply(amplitude_matrix.v,1,calcp,ratio_out$meanChangeRatio)
tpvalue = tpvalue[names(gsea_input)]
xnames = names(tpvalue)
trans_initialA = xnames[substr(xnames,1,3) == 'EIF']
trans_initialA = sort(trans_initialA)
trans_initialA
## [1] "EIF1AX" "EIF2A" "EIF2AK2" "EIF2B1" "EIF2B2" "EIF2B3"
## [7] "EIF2B4" "EIF2B5" "EIF2D" "EIF2S1" "EIF2S2" "EIF3A"
## [13] "EIF3B" "EIF3D" "EIF3E" "EIF3F" "EIF3G" "EIF3H"
## [19] "EIF3I" "EIF3J" "EIF3K" "EIF3L" "EIF3M" "EIF4A1"
## [25] "EIF4A2" "EIF4A3" "EIF4B" "EIF4E" "EIF4E2" "EIF4EBP1"
## [31] "EIF4G1" "EIF4G2" "EIF4G3" "EIF4H" "EIF5A" "EIF5A2"
## [37] "EIF5B" "EIF6"
trans_initialB = xnames[substr(xnames,1,3) == 'EEF']
trans_initialB = sort(trans_initialB)
trans_initialB
## [1] "EEF1A1" "EEF1A2" "EEF1AKMT1" "EEF1B2" "EEF1E1" "EEF1G"
## [7] "EEF2"
trans_initialA = c(trans_initialA,trans_initialB)
for(i in 1:length(pro_amplitude)){
tx = pro_amplitude[[i]][trans_initialA]
#tx= tx[!is.na(tx)]
tmp = data.frame(genes = trans_initialA,stringsAsFactors = F,
amplitude = as.vector(tx))
if(i ==1){
tmpout = tmp;
}else{
tmpout = cbind(tmpout,tmp[,'amplitude'])
}
}
tmpout = t(tmpout[,-1])
rownames(tmpout) = names(pro_amplitude)
colnames(tmpout) = trans_initialA
idx = colSums(is.na(tmpout)) < 3
sum(idx)
## [1] 31
#tmpout = tmpout[rowSums(is.na(tmpout)) < 3,]
tmpout = tmpout[,idx]
metadata <- tissueClass[rownames(tmpout),c('type','type')]
tclass = data.frame(class = metadata$type,row.names = rownames(metadata))
colnames(tclass) <- c("Tissue_type")
ann_colors = list(
Tissue_type= c('#FD60A7','gray85','#4fffF7')
)
names(ann_colors$Tissue_type) = c('Type I','Undefine', 'Type II')
pheatmap::pheatmap(t(tmpout),annotation_col = tclass,cluster_rows = F,
annotation_colors = ann_colors,
main = 'Change in initiation and elongation factors',
color=colorRampPalette(c('#3B4992','gray99','#EE0000'))(30),
#file ="./out/20230217_aging/Figure6_heatmap_TED_machanism/Figure6B_heatmap_translation_elongationV3.pdf",
height = 8,width = 8
)
11.4 Figure 7C
xnames = names(tpvalue)
ribosome = file2frame('./data/ribosome_proteins_from_kegg.txt',header = F)
ribosome = ribosome$V1
trans_initialA = intersect(ribosome,xnames)
trans_initialA
## [1] "RPS2" "RPS3" "RPS3A" "RPS4X" "RPS5" "RPS7" "RPS9" "RPS13"
## [9] "RPS14" "RPS15" "RPS16" "RPS20" "RPS21" "RPS23" "RPS24" "RPS27L"
## [17] "RPS27A" "RPS29" "FAU" "RPSA" "RPL4" "RPL5" "RPL6" "RPL7"
## [25] "RPL7A" "RPL8" "RPL12" "RPL13" "RPL13A" "RPL14" "RPL18" "RPL22"
## [33] "RPL23" "RPL23A" "RPL27" "RPL28" "RPL30" "RPL32" "RPL35" "RPL38"
## [41] "RPLP1" "RPLP2"
for(i in 1:length(pro_amplitude)){
tx = pro_amplitude[[i]][trans_initialA]
#tx= tx[!is.na(tx)]
tmp = data.frame(genes = trans_initialA,stringsAsFactors = F,
amplitude = as.vector(tx))
if(i ==1){
tmpout = tmp;
}else{
tmpout = cbind(tmpout,tmp[,'amplitude'])
}
}
tmpout = t(tmpout[,-1])
rownames(tmpout) = names(pro_amplitude)
colnames(tmpout) = trans_initialA
idx = colSums(is.na(tmpout)) < 3
sum(idx)
## [1] 34
#tmpout = tmpout[rowSums(is.na(tmpout)) < 5,]
tmpout = tmpout[,idx]
tmpout = tmpout[rownames(tissueClass),]
metadata <- tissueClass#[rownames(tmpout),]
#metadata <- tissueClass[rownames(tmpout),c('class','class')]
tclass = data.frame(class = metadata$type,row.names = rownames(metadata))
colnames(tclass) <- c("Tissue_type")
ann_colors = list(
Tissue_type= c('#FD60A7','gray85','#4fffF7')
)
names(ann_colors$Tissue_type) = c('Type I','Undefine', 'Type II')
pheatmap::pheatmap(t(tmpout),annotation_col = tclass,cluster_rows = F,
annotation_colors = ann_colors,
main = 'Change in ribosomal proteins',
cluster_cols = T,color=colorRampPalette(c('#3B4992','gray99','#EE0000FF'))(30),
#file ="./out/20230217_aging/Figure6_heatmap_TED_machanism/Figure6C_heatmap_ribosomeV3.pdf",
height = 8,width = 8
)
